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Abstract 

In this article, an on-line EM algorithm is derived for general Exponential Family models with 
Hidden variables (EFH models). It is proven that the on-line EM algorithm is equivalent to a 
stochastic gradient method with the inverse of the Fisher information matrix as a coefficient matrix. 
As a result, the stochastic approximation theory guarantees the convergence to a local maximum 
of the likelihood function. 

The performance of the on-line EM algorithm is examined by using the mixture of Gaussian 
model, which is a special type of the EFH model. The simulation results show that the on-line EM 
algorithm is much faster than the batch EM algorithm and the on-line gradient ascent algorithm. 
The fast learning speed is achieved by the systematic design of the learning rate schedule. Moreover, 
it is shown that the on-line EM algorithm can escape from a local maximum of the likelihood 
function in the early training phase, even when the batch EM algorithm is trapped to a local 
maximum solution. 
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1 Introduction 

In recent years, the EM algorithm (Dempster, Laird and Rubin 1977) has been applied to a number 
of neural network models (Jordan & Jacobs 1994: Amari 1995; Xu et al. 1995). The EM algorithm is 
a general method for finding a maximum likelihood estimator for a stochastic model with hidden vari¬ 
ables. Many neural network models have corresponding stochastic models (Amari 1995: Bishop 1995). 
The stochastic formulation helps to study the learning problem of neural networks. Hidden variables 
in the stochastic models correspond to hidden units in the neural networks, and these variables expand 
the representational capability of the models. 

The EM algorithm is a batch algorithm in which all the training data are used to update model 
parameters at each iteration. It was proven that the likelihood for the training data always increases 
(or does not change) at each iteration of the EM algorithm (Dempster et al. 1977). Therefore, the 
EM algorithm guarantees the convergence to a local maximum of the likelihood function. 

In real world problems, it is often desirable to use an on-line learning algorithm in which the 
training data are supplied one by one and the model parameters are updated each time using the 
current data. The most popular on-line algorithm is the on-line gradient descent algorithm. This 
algorithm can be considered as a stochastic approximation method (Robbins & Monro 1951) for 
finding a local minimum of an error function (Amari 1967; Kushner & Yin 1997; Bottou 1999; Murata 
1999). On-line EM algorithms have also been proposed by several authors (Jordan & Jacobs 1994 
; Nowlan 1991). However, there has been no theoretical study on the convergence of the on-line 
EM algorithms. Neal and Hinton (1998) proposed a wide variety of incremental EM algorithms and 
proved their convergence. However, one drawback of their study is that their algorithms either need 
additional storage variables for all of the training data or need to see all of the training data at each 
iteration. This prevents these algorithms from being applied to real-time on-line settings where new 
data are indefinitely supplied each time. 

In our previous paper, we proved the convergence of the on-line EM algorithm for a normalized 
Gaussian network (Sato & Ishii 1998). In this paper, we study more general models, namely the 
Exponential Family models with Hidden variables (EFH models), and derive an on-line EM algorithm 
for them. We prove that the on-line EM algorithm is equivalent to a stochastic gradient method 
with the inverse of the Fisher information matrix as a coefficient matrix. As a result, the stochastic 
approximation theory guarantees the convergence to a local maximum of the likelihood function. 

The performance of the on-line EM algorithm is examined by using the mixture of Gaussian model, 
which is a special type of the EFH model. The simulation results show that the on-line EM algorithm 
is much faster than the batch EM algorithm and the on-line gradient ascent algorithm. The fast 
learning speed is achieved by the systematic design of the learning rate schedule. Moreover, it is 
shown that the on-line EM algorithm can escape from a local maximum of the likelihood function in 
the early training phase, even when the batch EM algorithm is trapped to a local maximum solution. 

The paper is organized as follows. Section 2 explains three kinds of learning problems for stochas¬ 
tic models. In this paper, only unsupervised learning problem is studied because the extension of 
the current method to the two kinds of supervised learning problems is obvious. Section 3 is a short 
review of the EFH model. The on-line EM algorithm for the EFH model is derived in Section 4. The 
equivalence of the on-line EM algorithm to the stochastic gradient method is proven in Section 5. Sec¬ 
tion 6 explains the systematic design of the discount factor schedule. Section 7 explains experimental 
results and Section 8 is a conclusion of this paper. 


2 Learning Problem 

In this section, three kinds of learning problems are defined for stochastic models. 
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2.1 Unsupervised Learning 

For a set of observed input data X = {x(t)\t = 1. , unsupervised learning can be defined as a 

problem to find the stochastic model that best explains the observed input data. A stochastic model 
for the input distribution is defined by a probability distribution for an input x, P(x\6). where 6 
represents a set of model parameters. Unsupervised learning for the stochastic model can be solved 
by the maximum likelihood method, which maximizes the log-likelihood for a set of observed data 
X, logP(a;(i)|#). This model tries to learn the input data distribution p(x). The mixture of 

Gaussian model is an example of this type. 

2.2 Supervised Learning 

Deterministic neural networks are general purpose function approximators with adjustable (weight) 
parameters (Bishop 1995). For a set of observed input data X = {x(t)\t = 1 and the corre¬ 

sponding target data (teacher signal) Y = {y(t)\t = 1, ...,T}, supervised learning can be defined as a 
problem to find the optimal model parameter that minimizes some error function for a set of observed 
data (X,Y). 

Many neural network models have corresponding stochastic models (Amari 1995; Bishop 1995). 
There are two kinds of supervised learning for stochastic models. The first kind of supervised learning 
defines a stochastic model by a conditional probability for an output y given an input x , P(y\x,6). 
Supervised learning can be defined as a problem to find the optimal parameter that maximizes the 
log-likelihood of the conditional probability for a set of observed data (X, Y). l°g T > (y(t)|a;(f), 6). 
This model tries to learn the input-output relationship. In this case, the model does not take into 
account the input data distribution. Mixtures of expert models are examples of this type (Jacobs et 
al. 1991; Jordan & Jacobs 1994). 

After learning, a deterministic output for a given input x can be calculated as the expectation 
value of the output: 

V = f dn(y)yP{y\x,0), (1) 

where dp(y) denotes a measure on the output space. 

The second kind of supervised learning defines a model by a joint probability for an input x and 
an output y, P(x,y\6). Here, supervised learning can be defined as a problem to find the optimal 
parameter that maximizes the log-likelihood of the joint probability, YlJ=i log 2 /(t)|@). This 
model tries to learn the joint input-output distribution. Therefore, the model takes into account the 
input data distribution as well as the input-output relationship. The normalized Gaussian network 
(Sato & Ishii 1998) is an example of this type. The conditional probability of the model can be 
calculated as 

P(ty\x,6) = P(x,y\,0)/ J dp(y')P(x,y'\0). (2) 

The deterministic output can be calculated by using Eq. (1) and (2). 

One can define the Exponential Family models with Hidden variables (EFH models) for the above 
three types of stochastic models and derive an on-line EM algorithm for each type. In this paper, 
however, we only study the on-line EM algorithm for unsupervised learning. This is done purely for 
simplicity. It is obvious that the on-line EM algorithm for the two kinds of supervised learning can be 
derived in the same way as in this paper and that one can also prove the convergence of the on-line 
EM algorithm for these models. 
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3 Exponential Family model with Hidden variables 

Let us define an Exponential Family model with Hidden variables (EFH model) for the input data 
distribution (Amari 1985). The input variable x = {x \,..., xd) T could be a vector in the D-dimensional 
Euclidean space or could be a discrete variable vector. The hidden variable is assumed to be a discrete 
variable vector and is denoted by z = (zi ,..., zl) t . A complete event of the model is specified by 

(x, z). An observed event x is called an incomplete event (Dempster et al. 1977). An EFH model is 

defined by a probability distribution for an observed event x : 

P(x|0) = ^P(x,z|0), (3a) 

{z} 

M 

P(x.z|6>) = exp[^] D'( x , z)0j + rM+i (x, z) T(0)], (3b) 

i =l 

where 0 = {0 \..... ()\i ) 7 denotes a set of model parameters and J2{z} denotes a sum over possible 
configurations for z. P(x, z\0) in Eq. (3) represents the probability distribution for a complete event 
(x, z). A set of sufficient statistics is defined by 

r(x,z) = (ri(x,z),...,r M (x,z)) T . (4) 

The normalization factor H/(0) is determined by 

exp[^(0)] = [ d/r(x)^exp[r T (x, z)-0 + r Ai f + i(x, z)], (5) 

' M 

where d/r(x) is a measure on the input space and r T ■ 6 = X)f=i r i@i- Equation (5) is derived from a 
probability condition: 

f d/j,(x.) ^ P(x, z\0) = l. (6) 

{Z} 

The expectation parameter for the EFH model. 4> = {<j> i 5 (Amari 1985) is defined by 

</) = d^(O)/d0 = (d^/d6 1 ,... : d^/dd M ) T (7a) 

= E[r(x,z)\0] = I dfx(x) ^r(x,z)P(x,z|0). (7b) 

{z} 

It is known that the transformation between 0 and is one-to-one. Therefore, Eq. (7) can be inverted 
as 

0 = 0 ( 0 )- ( 8 ) 

The expectation parameter (f) can be used to parameterize the probability distributions of the EFH 
model and it is a dual set of the model parameter 0 (Amari 1985). By using the negative of the 
entropy, 

H{(j)) = I <:fyi(x)^P(x,z|0(0))logP(x,z|0(<£)) (9a) 

{z} 

= 4> t •0-'F(0), (9b) 

the model parameter 0 can be explicitly expressed as 


0 = e{<t>) = dH(4>)/d<f). 


( 10 ) 
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Equations (7), (9) and (10) form the Legendre transformation. 

The Mixture of Exponential Family (MEF) model is a special type of the EFH models. The MEF 
model consists of N units. It is assumed that a single unit is selected from a set of N units for each 
observed event x. The hidden variable of the MEF model is an indicator variable, z = (zi ,..., zn) t . 
If the n-th unit is selected, then z n = 1 and z m = 0 for m / n. There are N configurations for the 
indicator variable z corresponding to the N units. The probability distribution for a complete event 
(x, z) in the MEF model is given by 


N Ii N 

P(x, z|0) = exp[^ Y z n r n (x)fl«.„ + Y Z <finX) + r K + l(x) - ^(0)], (11) 

n I Q=1 n=1 

where 0 = {0 rt .„|n = l,....N:a = 0,1,..., K} is a set of MEF model parameters. The sufficient 
statistics of the MEF model are defined by r(x, z) = {z n r Q (x)\n = 1 ,....N:a = 0,where 
ro(x) = 1 is assumed. Using Eq. (3a) and (11), one can obtain the MEF model distribution for an 
observed event x: 

N 

P(x|0) = £P(x,n|0) (12a) 

71=1 

K 

P(x, n\6) = exp[J^ r a (x)6 n , a + 9 nfi + r K+ i(x) - ^(0)] {n = 1,..., N). (12b) 

a =1 


P(x,n|0) represents a joint probability for a complete event such that the n-th unit is selected and 
the input x is observed. 

For later convenience, let us define the Fisher information matrices G(0) and V($) for the prob¬ 
ability distributions P(x, z|0) and P(x.z|0(0)), respectively: 


Gij(O) 

ViM) 


E 


E 


/51ogP(x, z|0)\ /31ogP(x,z|0 


V 


do, 


) 


d0. 


/ aiogP(x.z|0(^)) ^ I fc)logP(x,z|0(^)) ^| 


d(j), 

(i.j = 1,..., M). 


'A 


d 4> j 


The following relations are important for later calculations: 


Gij(0) 


ViM) 

G (0) 


9 2 T(0) _ d(k 
dOjdOj dOj 
d 2 H(c/)) _ &Y 
d<f>id(/)j d<i>, 

v(0)- 1 . 


(13a) 

(13b) 

(13c) 

(14a) 

(14b) 

(14c) 


4 On-line EM algorithm 

From a set of observed data X = {x(t) | t = 1, ...,T}, the EFH model parameter 0 can be determined 
by the maximum likelihood method. In particular, the EM algorithm (Dempster, Laird and Rubin 
1977) can be applied to the EFH model. The EM algorithm repeats the following E-step and M-step. 
It was proven that the likelihood for a set of observations increases (or does not change) after an E- 
and M-step. Therefore, the EM algorithm guarantees the convergence to a local maximum or a saddle 
point of the likelihood function. 
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• E (Estimation) step 

Let 6 be the present estimator. By using 0, the posterior probability of the hidden variable z 
for each observation x(t) is calculated according to the Bayes rule. 

P(z|x(t),0) = P(x(t),z|0)/P(x(t)|0) (15a) 

P(x(f)|0) = ^2 P( x (t), z|0). (15b) 

{z} 

• M (Maximization) step 

By using the posterior probability (15). the expected log-likelihood Q(G\6.~X.) for the complete 
events is defined by 


T 

Q(0 |0,X) = ££P(z|x(i),0)logP(x(i),z|0). (16) 

t=l {Z} 

On the other hand, the log-likelihood of the observed data X is given by 

L(0\X) = £logP(x(t)|0) = £log (eWM*)] • (17) 

t=i t=i \{z} ) 

Since an increase of Q(6\0, X) implies an increase of the log-likelihood L(0 |X) (Dempster, Laird 
and Rubin 1977), Q(6\G,'X.) is maximized with respect to the estimator 6. The stationary 
condition dQjdO = 0 can be solved as 


(t> = (r(x,z))(T). (18) 

The symbol (•) denotes a weighted mean with respect to the posterior probability (15) and is 
defined by 

(r(x,z))(T) = ^££r(x(i),z)P(z|x(i),£). (19) 

t=i {Z } 

The new model parameter 0 can be obtained by using Eq. (10): 

0 = 0(0) = SP(0)/d0. (20) 

The above EM algorithm is based on batch learning, namely, the parameters are updated after 
seeing all of the observed data X. In the following, we derive an on-line version of the EM algorithm. 
In the case of on-line learning, the observed data {x(l),x(2),...} are supplied one at a time, and the 
estimator is changed after each observation. Let 6(t) be the estimator after the t -th observation x(f) 
and introduce a time-dependent discount factor A(i)(0 < A(f) < 1; t = 2, 3,...). A discounted weighted 
mean is defined by 


<C r(x, z) » ( t ) 


v(t) 


V(t) £ II A (' s ) 

T= 1 \S=T+1 

e ri a(.) 

lT=1S=T+1 


£ r ( x (T),z)P(z|x(T),0( T - 1)) 
{z} 

(t = 1,2,...)- 


(21a) 

(21b) 
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The discount factor A (t) is introduced for ignoring the effect of the old posterior values that employ 
the earlier inaccurate estimator. The schedule of the discount factor is discussed in Section 6. rj(t) is 
a normalization coefficient and plays the role of a learning rate. Replacing the weighted mean {•) by 
the discounted weighted mean <C • 3>. the on-line EM algorithm is obtained. 

The discounted weighted mean • >> can be calculated by a step-wise equation: 


< r(x. z) » (t) = «r(x,z)»(f —l)+?j(f) 


^r(x(t),z)P(z|f)- < r(x, z) > (t 

Lfz} 


r]{t) = (l + X(t)/r](t - 1)) \ 

where P(z|t) = P(z|x(t),0(t — 1)). Subsequently, the new parameters are calculated as 

4>{t) = < r(x, z) > (t) 

»(*) = «<*«> = ^rVr 

Equations (22) and (23) define the on-line EM algorithm for the EFH model. 


1) 


(22a) 


(22b) 


(23a) 

(23b) 


5 Stochastic approximation 


If an infinite number of data, which are drawn independently according to the input data distribution 
p(x), is available, the log-likelihood function is given by 

P(0) = P[log£P(x,z|0)]„ (24) 

fz} 

where E[-] p denotes the expectation value with respect to the input data distribution p(x). The 
maximum likelihood condition 8L(6)/d6 = 0 turns out to be 

4> = £[£r(x,z)P(z|x,0)] p (25a) 

M 

e = 0(0) = aff(0)/00. (25b) 

The on-line EM algorithm, Eq. (22) and (23), can be written as 

A 0(f) = 0(f)-0(f-l) 

= r ? (t)E r (x(t),z)P(z|x(t),0(t - 1)) - 0(i - 1)] 
fz} 


em = = !*,„• 

Using the log-likelihood for the current data x(t), 

P(x(t)|6>) = log(^P(x(i),z|0)), 


fz} 


the update rule (26) is further rewritten as 

' dL(x(t)\0 


A0(f) = v(t) ( 


86 


0{t i; 

(<90 ^ ( SP(x(t)|0(0))^ 

v{t) M -J 


= ^)(vw.- 1 ) ))- i ( at( y )) )i^ 1 „ 


(26a) 

(26b) 


(27) 

(28a) 

(28b) 

(28c) 
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where the relation (14) is used. 

Equation (28) shows that the on-line EM algorithm (26) is equivalent to a stochastic gradient 
method with the inverse of the Fisher information matrix V -1 as a coefficient matrix. In order to 
apply the stochastic approximation theory (Kushner and Yin 1997), it is assumed that the expectation 
parameter (f)(t) is bounded in a compact region \4>i{t)\ < </> max (i = 1 ,...,M) for a sufficiently large 
constant (fi ma x- This can be achieved by a simple procedure where is set to tymaxi—tymax)-, if 
4>i{t) becomes larger (smaller) than 4> ma x(~<t>max)- In numerical calculations, this kind of procedure 
is often necessary to avoid numerical explosion. It is also assumed that the input data distribution 
p(x) has a compact support or that p(x) decreases exponentially as |x| goes to infinity. This is also a 
natural assumption because actual data are always distributed in a finite domain. These assumptions 
guarantee that 

E[|V- 1 (0)(SL(x|0(0))/a < ^)| 2 ] p < 00 . (29) 

The effective learning rate rf(t)(> 0) is assumed to satisfy the condition 

OO OO 

r](t) = oo and yr/ 2 (t)< oo. (30) 

/ i t =l 

Under the above conditions, the stochastic approximation theory (Kushner and Yin 1997) guaran¬ 
tees that the asymptotic behavior of the estimator satisfying (28) is described by a continuous time 
differential equation 

and that (f)(t) converges to a limit point of the differential equation (31) or converges to a boundary 
point of the compact region. Equation (31) has the Lyapunov function (—L(d(<fi))), which always 
decreases its value over time, because the Fisher information matrix and its inverse are positive definite 
matrices. Therefore, the estimator converges to a local maximum of the likelihood function L(6 (</>)) 
or converges to a boundary point of the compact region. In addition to the above assumptions, if 
(V -1 (4>)(dL(0(<f>))/d(f>)) at the boundary of the compact region always points inward, the estimator 
(f>(t) converges to a local maximum of the likelihood function inside the compact region with probability 
one. 

Equation (28c) has a similar form as the natural gradient method proposed by Arnari (1998), 
which gives the optimal asymptotic convergence. However, there is an important difference. The 
Fisher information matrix in (28c) is defined for the probability distribution of a complete event 
P(x, z\0((/>)). On the other hand, the Fisher information matrix in the natural gradient method is 
defined for the probability distribution of an observed event P(x|0(</>)). It is beyond the scope of the 
current paper to study the relationship between the on-line EM algorithm and the natural gradient 
method. 

6 Schedule of Discount Factor 

The performance of the on-line EM algorithm depends on the schedule of the effective learning rate 
rj(t), which is controlled by the schedule of the discount factor A (t). In this section, the schedule of 
the discount factor A (t) is discussed in detail. It is convenient to introduce new parameters by 


t t 


n(t) 

e e n a ( s ) = w*) 

T=1S=T+1 

(32a) 

e(f) 

= l-A(f). 

(32b) 
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n(t) represents an effective number of averaged data that contributes to the discounted weighted mean 
<C r(x. z) (t) defined by Eq. (21). In particular, n(t) = t for A(i) = 1. If A (t) = \(const.) , then 
A J = (1 — ef ~ 0 for t > (1/e). Therefore, (1/e) corresponds to an effective length of the memory 
window. The contribution from the past data earlier than this window length can be neglected in the 
calculation of the discounted weighted mean r(x, z) ( t ). 

From the definition (32), n(t) takes the maximum and minimum values for A (t) = 1 and A (t) = 0, 
respectively. Therefore, the following relations hold. 


1 < n(t) < t 

1 > 77 (f) >1 /t. 

If this constraint is satisfied, Eq. (22b) can be solved for A (t) as 

X(t) = r,(t - l)(l/r,(t) - 1). 


(33a) 

(33b) 


(34) 


Therefore, one can determine the effective learning rate rj(t) without referring to the discount factor 
A(f) as long as the constraint (33) is satisfied. In practice, however, it is easier to control the discount 
factor \(t) than the effective learning rate 77 (f) as will be explained later. 

In order to see asymptotic behavior of the effective learning rate 77 (f), the recursive equation (22b) 
is rewritten as 

nft) — nft — 1) = 1 — eft)nft — 1). (35) 

If e(t) —> f _ ^(/3 > 1) as t —> cx), then (n(t) — nft — 1)) —> 1, is satisfied. This implies that n(t) —>• f+o(l), 
i.e., 

77 (f) 1 1// for t~P((3> 1). (36) 

If eft) —> 1 /( 7 1) and n(t) — > (Tt), then the relation, T = 1 — T /y, can be derived from (35) in the limit 

t —> 00 , i.e., 

for 6(0^1 /(it), (37) 

where T and 7 are related by 

0 < r = 7/(7 + 1 ) < 1 . (38) 

If r](t) 1 /(Ft) is satisfied, the stochastic approximation condition (30) is satisfied. 

In the experiments, we used a (1/t)- schedule for e(t): 


e{t) 


1 

(f - 2)7 + l/e 0 ’ 


(39) 


where e( 2 ) = eo is the initial value of eft) and 7 determines the decay rate of eft). If the initial value 
of r]ft) is given by 1 

7(!) = 7 o = e 0 (l + 7 ), (40) 

the effective learning rate rjft) also satisfies a ( 1 /t)- schedule, 


vft) 


1 

ft ~ l)r + l/r]o' 


r = 7/(7 + 1)- 


( 41 ) 


1 In the definition of <g r(x, z) 7> (<) and 77(f), (21), it is assumed that there is no initial value contribution for 
<C r(x, z) S> (f). In the following arguments, we assume that there is an initial value contribution for <C r(x, z) 7 > (t), 
i.e., r(x, z) 1§> ( 0 ) = </>( 0 ) and 77(1) = 770. This is equivalent to assuming that we start at t = to > 1 and <j>{to — 1 ) is 
regarded as the initial value. 
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On the other hand, if the initial value does not satisfy Eq. (40). r](t) does not satisfy the (1/t)- schedule 
(41). In the experiments, we used three independent constants, eo,? 7 o and 7 , to control the learning 
schedule. They have clear physical meaning. The initial condition 77 ( 1) = 70 implies that the initial 
value for the expectation parameter </>(0) can be regarded as an average of (1 /770 — 1) data points. 
Therefore, the value of 770 is determined according to the uncertainty of the initial estimator </>(0). 
The initial condition e(2) = eo implies that the effective length of the memory window in the early 
stage of learning is given by (1/eo). The decay rate 7 controls the asymptotic behavior of 77 (f) as in 
Eq. (37) and (38). Fast and stable learning can be obtained by a schedule with r]o ~ o(l), 0 < eo 1 
and 0 < 7 < 1. There are three phases which are characterized by different 77 (f) behaviors (Figure 
2). In the first phase, 1 < f < 1/eo, e(i) is nearly equal to eo and 77 (f) rapidly decreases to eo- In this 
period, a rough estimate of the parameters takes place very quickly. This period can be considered an 
initialization phase. If the initial value for </>(0) is a good estimator, this period is not necessary and 
it is possible to use a (1/t)- schedule for 77 (f), (41),in which 770 is given by 770 = eo(l + 7 ) < 1. In the 
second phase, 1/eo < f < 1 /(eoT )5 e(i) an( l VW are nearly equal to eo- This period is the main part 
of the learning process, in which a fairly good estimator is obtained. In the third phase, l/(eo7) < ij 
77 (f) behaves like l/(Tf) in Eq. (37). This period stabilizes the estimator value by the (1/t)- schedule 
of 77 (f). 

In stochastic approximation theory, a schedule 77 (f) = 1/f is often used. This corresponds to 
A(f) = 1 in our on-line EM algorithm, i.e., there is no discount effect. Therefore, the effects of 
the old posterior values employing the earlier inaccurate estimator do not diminish and the initial 
parameter value greatly affects the learning performance. Nevertheless, asymptotic convergence to a 
local maximum of the likelihood function is guaranteed because contributions from an infinite number 
of data in the later stage of learning overwhelm the earlier finite number of contributions employing 
an inaccurate estimator. 


7 Experiment 

7.1 Mixture of Gaussian Model 

The performance of the on-line EM algorithm was examined by using the Mixture of Gaussian (MG) 
model. An MG model for the 2-dimensional input vector x = (aq,^) is defined by 

N 

P(x|6>) = ^P(n)exp[-|x-m„| 2 /(2cr2)](27rcr^)“ 1 . (42) 

n=l 

The center and the variance of the n-th Gaussian function are denoted by m,j and <x 2 , respectively. A 
set of model parameters is given by 6 = {0 n ,o, ni n , cr 2 \n = 1,.... N}, where P(n) = e e "’°/(Em= 1 

For a data generation model, an MG model consisting of four Gaussian functions was prepared. 
The centers and the variances of the model are shown in Fig. 1A. The mixing constant was given 
by (P(l),..., P(4)) = (4/9, 2/9, 2/9,1/9). Three training data sets consisting of 100, 1,000 and 10,000 
points were randomly generated according to the data generation model distribution p(x). 100 and 
1,000 data points in the training data sets are shown in Figs. 1A and IB, respectively. In order 
to measure the generalization error, the Kullback-Leiber (KL) divergence from the data generation 
distribution p(x) to the MG model distribution P(x|0) was calculated: 

K(p\\P) = J dx\dx 2 p('x.) Iog(p(x)/P(x|0)). (43) 

This was evaluated on a 51 x 51 point mesh grid. 
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The MG model was trained by using three algorithms: the on-line EM, the batch EM and the 
on-line gradient ascent algorithms. The on-line gradient ascent algorithm is defined by 

(44) 

For each learning method and each dataset, learning processes starting from 20 different initial con¬ 
ditions were performed. The initial conditions were generated as follows. First, the initial center 
positions were randomly generated from the region [0 1] x [0 1]. The initial variances for all the units 
were set to the variance among the initial centers. The mixing constants were set to the same value, 
(1/N). The same set of initial conditions was used for all experiments. 

In one epoch, each training algorithm sees all training data once. The batch EM algorithm updates 
the model parameters once at the end of each epoch, while the on-line algorithms update them for 
each datum. The computational costs in one epoch for three algorithms are on the same order that is 
proportional to (N x (Size of training data) ). The batch EM algorithm can be considered a special 
type of the on-line EM algorithm, in which the model parameters are only updated at the end of 
each epoch, and the discount factor is given by A (t) = 0 at the beginning of each epoch, or A (t) = 1 
otherwise. This fact also supports the assumption that the learning speed of each method can be 
compared in terms of number of epochs. 

7.2 Schedule of Discount Factor 

In order to achieve fast and stable learning performance, we used a (1/t)- schedule for e(t), (39). In 
the current experiments, the initial parameter value was randomly chosen. Therefore, we set rjo = 0.5, 
which means that the initial parameter is regarded as a contribution from one data point. The initial 
value eo is set to 0 . 01 , which means that a weighted average of the last 100 data points are taken 
in the early stage of learning. The decay rate 7 is set to 0.05. The performance of the on-line EM 
algorithm in the current experiments are fairly robust for the variation of these constants if rjo rsj o(l), 
0.001 < eo < 0.01 and 0.01 < 7 < 0.1 are satisfied. 

The time course of r](t) for various schedules are plotted in Fig. 2. A schedule r\(t) = 1/t, which 
corresponds to r/o = 1 and e(t) = 0, exhibits a fast decrease of r](t). The learning curves for this 
schedule are similar to those for the batch EM algorithm shown in Fig. 4. In a (1/t)- schedule for 
rj(t), (41), with 7 ]0 = 0.5 and 7 = 0.1 (eo = 7 o/(l + 7 ) — 0.45), r](t) decreases very slowly in the early 
stage of learning. The learning performance in this schedule is not good because the estimator forgets 
the contribution of past data too fast. The (1/t)- schedules for e(f), (39), with 70 = 0.5, eo = 0.01, and 
7 = {0.1,0.05,0.01} exhibit a fast decrease of r](t) in the early stage of learning and a slow decrease 
of rj(t) in the later stage. These schedules worked well in the experiments. 

If the initial parameter value is a good estimator, one can use a (1/t)- schedule for rj(t) with a 
small initial value of 70 = eo (1 + 7 ) — eo 1 . 

7.3 Results 

The learning curves for the on-line EM algorithm are shown in Fig. 3. The estimator quickly converged 
close to the optimal value in all cases. It nearly converged before 20,000 data points were presented, 
regardless of the training data size. The final generalization error improved as the amount of training 
data increased. Although over-training can be seen for the data set with 100 data (Fig. 3A), the 
generalization error at the asymptotic regime is the same as that obtained by the batch EM algorithm 
(see later discussion). The MG model with 20 units has an excessive degree of freedom. Therefore, 
there are a continuum set of optimal estimators and a continuum set of suboptimal estimators. This 
explains the small variations in the generalization error for the case N=20 (Fig. 3D). 
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The learning curves for the batch EM algorithm are shown in Fig. 4. Surprisingly, the learning 
speed of the batch EM algorithm was much slower than that of the on-line EM algorithm. The final 
generalization errors improved as the amount of training data increased and they were almost the 
same as those obtained by the on-line EM algorithm, except in the exceptional cases described below. 
The learning speed, however, did not improve even if the amount of training data increased. The 
estimator of an MG model with N=4 was trapped to a local maximum of the likelihood function for 
some initial conditions (Figs. 4A-4C). The increase in the training data size did not help this situation. 
The configuration of this local maximum solution is shown in Fig. 1C. In this solution, two units try 
to approximate one Gaussian function, so that the remaining two units must try to approximate 
three Gaussian functions. The estimator trained by the on-line EM algorithm nearly converged to the 
global maximum likelihood estimator, even when the same initial condition was used. The solution 
obtained by the on-line EM algorithm is shown in Fig. ID. If an MG model with 20 units was used, 
the estimator seemed to converge to suboptimal estimators (Fig. 4D), as in the case of the on-line 
EM algorithm. 

The reason for the slow learning and the local maximum problem can be interpreted as follows. 
In the early stage of the batch EM learning, the same posterior probability employing an inaccurate 
estimator is used for all training data in one epoch. This may cause an inappropriate assignment of 
the units for the training data and result in a poor estimation of the new parameter. In the on-line 
EM algorithm, the posterior probability is calculated with the new parameter, which is updated each 
time data is presented. Therefore, the unit assignment gradually improves as more data are presented. 

The learning curves for the on-line gradient ascent algorithm are shown in Fig. 5, in which the 
best learning rate schedule is used. The learning performance of the on-line gradient ascent algorithm 
was very sensitive to the learning rate schedule. If the initial value for the learning rate rg, was not 
smaller than 0.001, unstable oscillatory behavior occurred in the early stage of learning. The best 
performance was obtained by a (1/t)- schedule for rj(t), (41), with r/o = 0.0001 and 7 = 0.1, (i.e., 
eo = 7 o/(l + 7) — 0.00009, F = 7/(1 + 7) ~ 0.09). The performance for the fixed learning rate, 
r](t) = r/o = 0.0001, was slightly worse than the above best performance. In order to achieve the same 
accuracy of the generalization error obtained by the on-line EM algorithm, 100 times more epochs are 
needed in the on-line gradient ascent algorithm for all cases in Fig. 5. 

8 Conclusion 

In this article, an on-line EM algorithm was derived for general Exponential Family models with 
Hidden variables (EFH models). It was proven that the on-line EM algorithm is equivalent to a 
stochastic gradient method with the inverse of the Fisher information matrix as a coefficient matrix. 
As a result, the stochastic approximation theory guarantees the convergence to a local maximum of 
the likelihood function. 

The performance of the on-line EM algorithm was examined by using the mixture of Gaussian 
model. The simulation results showed that the on-line EM algorithm was much faster than the batch 
EM algorithm and the on-line gradient ascent algorithm. The efficiency of the on-line EM algorithm 
over the batch EM algorithm became more prominent as the amount of training data increased. 
In addition, the on-line EM algorithm could escape from a local maximum, even when the batch 
EM algorithm was trapped to a local maximum solution. The inherent stochastic nature of the on¬ 
line EM algorithm in the early stage of learning may be the reason for this fact. This also shows 
the superiority of the on-line EM algorithm over the batch EM algorithm, although this does not 
guarantee the convergence to the global maximum. 

One of the advantages of our on-line EM algorithm is that the estimator at each time step can 
be represented by the discounted weighted mean <C r(x, z) (t) defined by Eq. (21). Therefore, 
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the on-line EM algorithm gives a reasonable estimator even in the early stage of learning. Moreover, 
the learning rate schedule can be systematically designed and the fast learning schedule, in which 
r]Q ~ o(l), is possible. This kind of fast learning schedule is not possible for the on-line gradient ascent 
algorithm because it causes unstable oscillatory learning behaviors. 

We have pointed out that the on-line EM algorithm has a similar form as the natural gradient 
method proposed by Amari (1998), which gives the optimal asymptotic convergence. The inverse of the 
Fisher information matrix in the on-line EM algorithm may contribute to fast learning performance. In 
our on-line EM algorithm, however, it is not necessary to calculate the inverse of the Fisher information 
matrix. In the future, it would be interesting to study the relation of our algorithm to the natural 
gradient method. 
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• Figure 1 

(A), (B): Centers and variances of the data generation model consisting of four Gaussian func¬ 
tions are shown by four circles. The center and radius of a circle correspond to the center and 
standard deviation of a Gaussian function, respectively. 100 and 1,000 training data points are 
plotted in (A) and (B), respectively. (C): Local maximum solution obtained by the batch EM 
algorithm is shown by four circles, which correspond to four Gaussian functions in the trained 
MG model. (D): Solution near the global maximum likelihood estimator obtained by the on-line 
EM algorithm is shown by four circles. The same initial condition is used in (C) and (D). 

• Figure 2 

The time courses of the effective learning rate 77 (f) for various schedules are plotted. Dashed 
and dash-dotted lines represent a schedule 77 (f) = 1/f and a (1/t)- schedule for 77 (f), (41), with 
770 = 0.5, and 7 = 0.1, respectively. Solid lines represent the (1/t)- schedules for e(f), (39), with 
770 = 0.5, eo = 0.01, and 7 = {0.1,0.05,0.01}. 

• Figure 3 

Learning curves for the on-line EM algorithm. Ordinate denotes the generalization error mea¬ 
sured by the KL-divergence from the data generation model distribution to the trained MG 
model probability distribution. Abscissa denotes the epoch number. Each figure plots learning 
curves starting from 20 different initial conditions. The sizes of the training data are 100 (A); 
1,000 (B); 10,000 (C); and 1,000 (D). The numbers of units in the trained MG models are four 
(A, B, and C) and 20 (D). 

• Figure 4 

Learning curves for the batch EM algorithm. The training data sets and the trained MG models 
are the same as in Figure 3. 

• Figure 5 

Learning curves for the on-line gradient ascent algorithm. The training data sets and the trained 
MG models are the same as in Figure 3. 








ETA 


FIGURES 


17 



Figure 2 





FIGURES 



100 0 



50 

Epoch 


100 0 


Figure 3 






FIGURES 



100 0 



50 

Epoch 


100 0 


Figure 4 










